clear all
clf
tic

F1='201201-202203.xlsx';%�ļ����ַ��� 
A1=xlsread(F1);

A_in=A1(1:108,2)./1e8;
A_out=A1(1:108,3)./1e8;

ABTP=A_in;

L=60;
Tp=12;
BTP=ABTP(end-L+1:end,:);
M=48;
N=12;

% in alp=0.4;
% out alp=0.1;
alp=0.4;

Btp=BTP(1:M,:);
test_trabtp=BTP(M:M+N-1,:);
test_btp=BTP(M+1:M+N,:);

Ut=max(Btp)-min(Btp);

T=30;
mRMSE=zeros(1,T);
C_n_star=2;

for C_n=C_n_star:T

aU=Ut/C_n;

Utl=zeros(C_n,1);
Uth=zeros(C_n,1);

for i=1:C_n
    Utl(i,1)=min(Btp)+(i-1)*aU;
    Uth(i,1)=min(Btp)+i*aU;
end

center=(Utl+Uth)./2;
Au=zeros(1,M);

for i=1:M
    dist = distfcm(center, Btp(i,:));
    [~, Au(i)] = min(dist);% 
end

Au1=Au(1,1:M-1);
Au2=Au(1,2:M);

R=zeros(C_n,C_n);
for i=1:M-1
    R(Au1(i),Au2(i))=R(Au1(i),Au2(i))+1;
end

%R(find(R>0)) = 1;
w=R./sum(R,2);
w(find(isnan(w)==1)) = 0;

% ���ģ��ֵ
for i=1:C_n
    eval(['BAu_',num2str(i),'=find(Au==i);']);
    eval(['L_BAu_',num2str(i),'=length(BAu_',num2str(i),');']);  
    eval(['Sum_BAu_',num2str(i),'=sum(Btp(BAu_',num2str(i),',:));']);
end
%cen=center;
cen=zeros(C_n,1);
for i=1:C_n
    if i==1
       eval(['cen(i,1)=(Sum_BAu_',num2str(i),'+0.5*Sum_BAu_',num2str(i+1),')/(L_BAu_',num2str(i),'+0.5*L_BAu_',num2str(i+1),');']);
    elseif i==C_n
       eval(['cen(i,1)=(Sum_BAu_',num2str(i),'+0.5*Sum_BAu_',num2str(i-1),')/(L_BAu_',num2str(i),'+0.5*L_BAu_',num2str(i-1),');']);
    else
       eval(['cen(i,1)=(0.5*Sum_BAu_',num2str(i-1),'+Sum_BAu_',num2str(i),'+0.5*Sum_BAu_',num2str(i-1),')/(0.5*L_BAu_',num2str(i-1),'+L_BAu_',num2str(i),'+0.5*L_BAu_',num2str(i-1),');']);
    end    
end

cen(find(isnan(cen)==1)) = center(find(isnan(cen)==1));

F=sum(w.*cen,2);
F(find(F==0)) = cen(find(F==0));

Ft=F(Au);

Ft=alp.*Ft+(1-alp).*ABTP(end-L+1-Tp:end-L+1-Tp+M-1,:);

t_er=Ft-Btp;
mRMSE(C_n)=sqrt(sum(t_er.^2)/size(t_er,1));
end

[index,aindex]=min(mRMSE(C_n_star:T));
aindex=aindex+C_n_star-1;
disp(['Index=',num2str(aindex)]);

%% ��֤ 

C_n=aindex;
aU=Ut/C_n;

Utl=zeros(C_n,1);
Uth=zeros(C_n,1);

for i=1:C_n
    Utl(i,1)=min(Btp)+(i-1)*aU;
    Uth(i,1)=min(Btp)+i*aU;
end

center=(Utl+Uth)./2;
Au=zeros(1,M);

for i=1:M
    dist = distfcm(center, Btp(i,:));
    [~, Au(i)] = min(dist);% 
end

Au1=Au(1,1:M-1);
Au2=Au(1,2:M);

R=zeros(C_n,C_n);
for i=1:M-1
    R(Au1(i),Au2(i))=R(Au1(i),Au2(i))+1;
end

%R(find(R>0)) = 1;
w=R./sum(R,2);
w(find(isnan(w)==1)) = 0;

% ���ģ��ֵ
for i=1:C_n
    eval(['BAu_',num2str(i),'=find(Au==i);']);
    eval(['L_BAu_',num2str(i),'=length(BAu_',num2str(i),');']);  
    eval(['Sum_BAu_',num2str(i),'=sum(Btp(BAu_',num2str(i),',:));']);
end
%cen=center;
cen=zeros(C_n,1);
for i=1:C_n
    if i==1
       eval(['cen(i,1)=(Sum_BAu_',num2str(i),'+0.5*Sum_BAu_',num2str(i+1),')/(L_BAu_',num2str(i),'+0.5*L_BAu_',num2str(i+1),');']);
    elseif i==C_n
       eval(['cen(i,1)=(Sum_BAu_',num2str(i),'+0.5*Sum_BAu_',num2str(i-1),')/(L_BAu_',num2str(i),'+0.5*L_BAu_',num2str(i-1),');']);
    else
       eval(['cen(i,1)=(0.5*Sum_BAu_',num2str(i-1),'+Sum_BAu_',num2str(i),'+0.5*Sum_BAu_',num2str(i-1),')/(0.5*L_BAu_',num2str(i-1),'+L_BAu_',num2str(i),'+0.5*L_BAu_',num2str(i-1),');']);
    end    
end

cen(find(isnan(cen)==1)) = center(find(isnan(cen)==1));

F=sum(w.*cen,2);
F(find(F==0)) = cen(find(F==0));

Sc=zeros(N,1);
for i=1:N
    dist = distfcm(center, test_trabtp(i,:)); %���������ģ��C-��ֵ����distfcm(center,data),data��center��ŷ�Ͼ���
    [~, Au_new] = min(dist);
    Sc(i,:)=F(Au_new);
end

Sc=alp.*Sc+(1-alp).*ABTP(end-L+1-Tp+M:end-L+1-Tp+M+N-1,:);

er=Sc-test_btp;
RMSE=sqrt(sum(er.^2)/size(er,1));
MAXE=max(abs(er));
ARE=sum(abs(er)./test_btp)/size(er,1);
disp(['RMSE=',num2str(RMSE)]);
disp(['MAXE=',num2str(MAXE)]);
disp(['ARE=',num2str(ARE)]);

figure(1)
subplot(3,1,1)
plot(C_n_star:T,mRMSE(C_n_star:T),'-bo');

subplot(3,1,2)
plot(test_btp,'-bo');
hold on
plot(Sc,'-rs');
legend('train','test');

subplot(3,1,3)
plot(test_btp-Sc,'-bo');

toc